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TECHNICAL PAPER 


ANALYTICAL AND NUMERICAL STUDIES OF THE THERMOCAPILLARY FLOW 
IN A UNIFORMLY ROTATING FLOATING ZONE 

I. INTRODUCTION 


The accessibility of the microgravity environment of space by NASA’s Space Shuttle has stimulated 
much research into the processing of materials in low gravity. Processing in such an environment holds poten- 
tial for the development of higher quality products than can be produced in ground-based laboratories, and 
for the development of completely new materials. Modem solid-state electronics is based on the properties 
of high quality crystals and there is substantial evidence for believing that superior crystals can be grown 
from a melt in the absence of flow in the melt. Reviews describing the effects of convection for several 
crystal growing systems have been given by Pimputkar and Ostrach [1] and Ostrach [2] . This paper is con- 
cerned with the fluid dynamical behavior of the floating zone crystal growing system in low gravity. 

The floating zone system for unidirectional single-crystal growth from a low-quality polycrystalline 
material is shown schematically in Figure 1. The ring heater maintains a liquid floating zone that bridges 
the gap between the two cylindrical solid crystalline rods. Surface tension at the free melt surface holds the 
liquid together. The polycrystalline rod is fed into the melt zone where it is melted while the other rod is 
withdrawn from the melt as a single crystal. Many materials of interest absorb impurities from their con- 
tainer walls and, hence, the great advantage of the floating zone system is that contamination can be reduced 
or even eliminated. 

AXIS 


I 



Figure 1. A schematic drawing of the floating zone system for unidirectional single 
crystal growth. Typical shape distortions are shown for the crystal-melt interfaces 
and for the free surface in the presence of gravity. 



Although crystals of pure substances like silicon are often grown using the floating zone system, 
there are also many applications of the system where another substance (a dopant) is added to produce a 
crystal with certain desired properties. Ideally, the dopant should be uniformly distributed throughout the 
crystal but this is not what is achieved. Variations of dopant concentration in both the radial and axial 
directions are invariably present [3] . For semiconductors these concentration variations lead to significant 
variations in electrical conductivity throughout the crystal which substantially degrade the electronic pro- 
perties for many applications. 

The major factors affecting the dopant distribution in the crystal are the temperature and flow fields 
in the melt. In general, the crystal solidification interface is not flat, and, as a consequence of the differential 
segregation occurring upon solidification, gravitational buoyancy produces radial concentration variations 
of the dopant at the crystal interface. Also, radial flows near the crystal interface produce radial dopant 
variations. Upon solidification these dopant variations lead to dopant variations in the crystal. Time depend- 
ent flows lead to striations of dopant concentration variation in the crystal. 

Operating a floating zone system in space eliminates buoyant convection, but gravity is not the only 
driving force for flow in the system. Since surface tension is a function of temperature, the axial tempera- 
ture distribution produced by the ring heater at the free surface sidewall of the melt drives a thermocapillary 
(or Marangoni) flow. Thermocapillary flow in a floating zone in the absence of rotation has been the subject 
of several experimental studies [4,5] (see also Fig. 8 of this paper) and theoretical (numerical) studies [6]. 
In general, the results show torroidal cells which fully penetrate the melt zone. 

The analytical work presented in this paper is an extension of the work by Smith [7] and Smith and 
Greenspan [8]. These workers examined the idea that the thermocapillary flow could be confined to a thin 
boundary layer close to the sidewall by means of a uniform rotation of the floating zone about its axis of 
symmetry. Such a rotation is achieved by rotating the crystals in the same direction at the same rate. Hence, 
in a rotating system in space there may be effectively no flow in the interior of the melt zone away from the 
sidewall. Smith [7] and Smith and Greenspan [8] considered a half zone, cold at one end and hot at the 
other, and tackled analytically the linearized thermocapillary problem in zero gravity and in the absence of 
crystal growth. They found that the weak, linearized thermocapillary flow is confined to a sidewall boundary 
layer. To extend this work to more realistic conditions, the present authors considered a full zone, cold at 
both ends and hot in the interior, and tackled analytically and numerically the linear thermocapillary prob- 
lem and numerically the nonlinear thermocapillary problem. 

In Section II, the model is described, the governing differential equations are presented, and the 
values chosen for the floating zone dimensions and the other parameters are given. In Section III, linearized, 
analytical solutions are found using singular perturbation theory and the various sidewall boundary layers 
described by Greenspan [9] for rotating fluids. In Section IV linear and nonlinear numerical solutions are 
presented. In Section V the analytical and linearized numerical results are compared, and in Section VI the 
linear and nonlinear flows are discussed. The conclusions and suggestions for future work are stated in 
Section VII. Mathematical details of the analytical work and a brief description of the numerical code are 
given in the Appendices. 

Mention should be made of the numerical work on the effects of rotation in a floating zone by 
Kobayashi and Wilcox [10], These workers were concerned with the flows driven by differential rotation 
of the solid endwalls. Thermocapillary driving was not included and hence their work is not relevant to the 
present paper. Differential rotation of floating zone systems has often been used to mix the melt and to 
reduce the effects of azimuthal asymmetries of the ring heater. The latter advantage is also present with a 
uniformly rotating system. 
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II. FORMULATION 


2.1 The Model and Governing Equations 

The floating zone geometry is taken to be a right circular cylinder of fluid as shown in Figure 2. 
The assumptions required for this to be true are discussed in Subsection 2.3. Cylindrical coordinates (r, 6, z) 
with the z-axis along the axis of symmetry of the fluid cylinder and the origin at the center of the cylinder 
are used. The top and bottom endwalls of the fluid cylinders are rigid boundaries located at z = ±L. The free 
surface sidewall of the cylinder is located at r = a and is assumed to be non-deformable (see Subsection 2.3). 
The endwalls are rotated uniformly with angular velocity, £2, about the axis of symmetry and are maintained 
at the melting point temperature, T m . The sidewall is maintained at t = t m + AT cos 7rz/2L. This tempera- 
ture boundary condition is the driving force for the flow. Gravity and the centrifugal force due to rotation 
are ignored (see Subsection 2.3). 

AXIS 


FREE LIQUID 
SURFACE 


Figure 2. The simplified floating zone configuration chosen for the present studies. 

The melt is a right circular cylinder of fluid with rigid plane endwalls and a 
free cylindrical sidewall. Cylindrical coordinates aligned with the axis of 
symmetry and with the origin at mid-depth in the melt are used. 

The velocity in the melt is denoted by v(u, v, w), the pressure by p, and the temperature by T. 
The fluid density is denoted by p, the kinematic viscosity by v, the thermal diffusivity by k, and the surface 
tension by a. The surface tension is taken to vary linearly with temperature a = a Q - y (T - T m ), where y is 

the temperature coefficient and o Q is the surface tension at the reference temperature, T m . 

The equations are non-dimensionalized using the half depth, L, the density at the reference tempera- 
ture, p Q , an as yet undetermined velocity, U, and the temperature difference, AT. The pressure is non- 

dimensionalized using p 0 UJ2L. The resulting dimensionless governing equations for steady, axisymmetric 
flow in a frame of reference rotating with the zone are: 
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( 1 ) 


e (uu r + wu z -^J= -P r +2v + E(u rr + -± + u zz - , 

e ( uv r + wv z + 7^) = " 2u + E ( v rr + “ + v zz > 

e(uw r + ww z ) = -p z + E ^ Wrr + — + w zz ^ , 

Pre(uT r + wT z ) = E (j rr + j- + T zz ) , 


( 2 ) 


(3) 


(4) 


u 


u r + 7 + w z° > 


(5) 


where a suffix denotes differentiation with respect to that quantity. The dimensionless parameters, e, E and 
Pr are defined below after equation (13). 

The boundary conditions on the endwalls are for no flow and a constant temperature. On the side- 
wall, there is no radial flow and no azimuthal stress, and the axial temperature distribution leads to an axial 
thermocapillary stress. In dimensionless terms: 


u = v = w = 0 


T = 0 


on z = ±1 


u = 0 , 


v r — = 0 , 


w f = - MT Z , 


on r = A 


irz 

T = cos — . 

2 


(6a,b,c) 


(7) 

( 8 ) 
(9) 

( 10 ) 

( 11 ) 


For an axisymmetric problem which is also symmetrical about mid-depth, the following must also be satis- 
fied: 


u = v = w r = Tj. = 0 , on r = 0 , 


( 12 ) 
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u z = v z = w = T z = 0 , on z = 0 . 


(13) 


The dimensionless parameters are defined as follows: the aspect ratio A = a/L, the Ekman number 

E = v/fiL 2 , the Rossby number e = U/fiL, the Prandtl number Pr = v/k, the Marangoni number M = 7AT//L1U, 
where p is the dynamic viscosity. 

A stream function is introduced defined as, 

1 

w = - — (ri//) r , u=\p z . (14) 

The global mass conservation condition is, 


\p = 0 , on r = A and z = ±1 


(15) 


2.2 Parameter Values 

In order to make these theoretical studies as relevant as possible to real floating zone systems, silicon 
was chosen for the melt material and realistic values were chosen for the zone dimensions. Table 1 lists the 
values of the physical properties of silicon used in this paper. These values are the same as those used by 
Smith [7] and Smith and Greenspan [8] and they are taken from Reference 3. Table 2 lists the values of the 
zone dimensions and other quantities. The value of 5°C for AT may appear unrealistically small when the 
high melting point of silicon (1410°C) is considered but this is not necessarily so; this point is discussed 
further in Section VII. Table 3 lists the values of the dimensionless parameters. Since, at this stage in the 
paper U has not been determined, the values of e and M cannot yet be calculated. However, for completeness 
the values calculated later (Section III) are included in Tables 2 and 3. The values of other dimensionless 
parameters to be discussed in Subsection 2.3 are also included in Table 3. 


2.3 The Cylindrical Geometry Assumptions 

The melt zone geometry is simplified to that of a right circular cylinder of fluid in uniform rotation 
(Fig. 2). The top and bottom endwalls were assumed to be plane. In a real floating zone system (Fig. 1), 
these boundaries are determined from a fuller set of equations and boundary conditions than are used here, 
and will, in general, be curved. However, the deviations from planarity can be made small in many cases. 

In general, a liquid wets its solid phase and hence the contact angle (measured from the solid through 
the liquid to the face surface) is acute. This is true for silicon. To satisfy this boundary condition on plane 
parallel endwalls together with the interior zero-gravity condition of constant pressure, an hour-glass con- 
figuration, rather than a cylinder is required. However, since the melt extends out to the radius of the solid 
crystal rods, the contact angle will be formed on the curve at this edge and the surface at contact can be 
parallel to the symmetry axis. In this case, the solution for the free surface is a cylinder. With plane parallel 
endwalls in the presence of gravity (gravity parallel to the axis of symmetry), a quantitatively different 
hour-glass shape will form to satisfy the wetting boundary condition and the hydrostatic pressure gradient 
(Fig. 1). 
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TABLE 1. PHYSICAL PROPERTIES OF LIQUID SILICON AT 1410°C 


Density 

p o 

2.5 gm/cm 3 

Dynamic Viscosity 

u 

0.88 x 10“^ gm/cm sec 

Kinematic Viscosity 

V 

0.35 x 10“^ cm^/sec 

Thermal Conductivity 

X 

0.32 Joule/cm sec C 

Thermal Diffusivity 

< 

0.15 cm 2 /sec 

Surface Tension 

0 

o 

720 dyne/cm 

Temperature Coefficient 
of Surface Tension 

y 

0.43 dyne/cm C 


TABLE 2. DIMENSIONS OF THE FLOATING ZONE AND OTHER EXPERIMENTAL PARAMETERS 


Depth 

2L 

2 cm 

Radius 

a 

1 cm 

Rotation Rate 

n 

1 rad/sec 

Temperature Difference 

AT 

5 C° 

Characteristic Flow 
Velocity 

U = irE 1 / 3 Y AT/2 y 

58 cm/sec 


TABLE 3. DIMENSIONLESS PARAMETERS 


Aspect Ratio, A 

a/L 

1 

Ekman Number, E 

v/QL 2 

0.35 x 10 -2 

Prandtl Number, Pr 

v / k 

0.023 

Marangoni Number, M 

YAT/pU 

4.2 

Rossby Number, e 

U/QL 

58 

Capillary Number, C 

pU/Oo 

4.5 x 10 -4 

Rotational Bond Number, B 

P {] 2 a 3 /a 

o ' o 

3.5 x 10 -3 
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Additional forces which can distort the cylindrical sidewall solution discussed above are the normal 
stresses due to viscosity and the centrifugal force. For the free surface to remain cylindrical, these stresses 
must be less than that due to surface tension. The ratios of the forces are measured by the capillary number 

C = /iU/a 0 and the rotational Bond number B = p Q f2^a^/a 0 , respectively, and for a flat cylindrical sidewall 
they should both be small. The values given in Tables 1 and 2 are used to derive, C = 4.5 x 10'^ and B = 
3.5 x 10 . Thus, for silicon and the dimensions chosen, a right circular cylinder of fluid is a good assumption. 


III. ANALYTICAL SOLUTIONS 


3.1 Linearization 

Analytical solutions were found by linearizing equations (1) through (4). The momentum and tem- 
perature advection terms are dropped. This implies that the Rossby number, e, is small. 


3.2 The Temperature Solution 

The temperature problem reduces to the solution of the two-dimensional conduction equation with 
the boundary conditions given by equations (7) and (11). The exact solution in non dimensional form is, 


I 0 (*r/2) 

I o 0rA/2) 


cos(irz/2) 


(16) 


where I Q denotes the modified Bessel function of zero order. Since, as a consequence of the linearization, the 

temperature equation is decoupled from the flow field, equation (16) is the complete temperature solution 
for all of the analytical studies given in this paper. 


3.3 The Interior Solution 

For this paper’s chosen conditions (Table 3), the Ekman number, E, is small and, hence, from equa- 
tions (1) through (3) one expects viscous boundary layers on the zone boundaries and an interior flow in 
which viscous effects are small. The authors proceed using singular perturbation theory and expand to 
successive orders in E. 

To leading order in E, equations (1), (2), (3), and (5) for the interior give: 


3pj 


_ + 2Vt = Ur = 

3r I I 9z 


9pj 9wj 


9z 


= 0 


(17a,b,c,d) 


Since there is no mechanism in this problem for establishing a radial pressure gradient, equations (17) reduce 
to: 
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Uj = Vj = Wj = pj = = 0 . 


(18a,b,c,d,e) 


Note that this interior solution satisfies all the endwall boundary conditions [equations (6)] and two side- 
wall boundary conditions [equations (8) and (9)] , but it does not satisfy the thermocapillary stress condi- 
tion [equation (10)]. 


3.4 The Leading Order E^ 3 Sidewall Boundary Layer 

The flow conditions close to the sidewall are now examined. Smith [7] and Smith and Greenspan 
[8] in their study of the half zone found that a sidewall boundary layer of nondimensional thickness E^/ 3 

could satisfy the sidewall boundary conditions. For this problem an E^ 3 layer is also considered. To formu- 
late the boundary layer equations and boundary conditions, a transformation to a dimensionless boundary 

layer coordinate, £ = E' 1 / 3 (A-r), is made. The governing equations for this layer are well known [9] and are 
given in Appendix A. Wj/3 ~ 0(1) is defined, where the suffix, 1/3, denotes an E 1 / 3 layer quantity. The 
other amplitude scales are then found to be: 

u 1/3 ~ E 1 / 3 , v , /3 ~ 0(1) , p 1/3 ~E 1 /3 f ^j/3 ~ E 1 / 3 . 

The E*/ 3 layer has Ekman layers [9] associated with it on the endwalls. The Ekman compatibility 
condition [9] when transformed and scaled yields, 

-w 1 1/r 3 v 1/3 

w 1/3 = ± 2 E If" ’ on z = ±1 ’ 


where the tilda denotes an 0(1) quantity. To leading order, therefore, wj^ is zero. Thus, the axial flow in 

the E^ 3 layer recirculates and does not enter the Ekman layer. The boundary condition for the E 1 / 3 layer 
becomes wj^ or ^1/3 = 0, on z = ±1. 

The axial thermocapillary stress boundary condition, equation (10), when transformed to the 
boundary layer coordinate and scaled becomes 


dwj/s 

~w~ 


— E^/ 3 M7rsin —nz 
2 2 


E*/ 3 try AT 1 

— — — sin — 7rz 

2pU 2 


Since the boundary layer thickness has now been determined, the characteristic velocity, U, can be deter- 
mined by balancing viscous shear against the thermocapillary stress. The equation 


U = 


E 1 / 3 iryAT 

2/i 
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is obtained and the thermocapillary stress boundary condition reduces to 


3 w 1/3 

a* 



= o 


on £ = 0 . 


The azimuthal stress boundary condition, equation (9), when transformed and scaled becomes. 


9 v l/3 

3 £ 


= 0 


on £ = 0 . 


When solving for the E ^ 2 layer, the problem is formulated in terms of \p j j 3 , and 


95 *1/3 

9£ 5 


on £ = 0 


is used for the azimuthal stress condition. This condition is equivalent to 


9 9 v 1/3 

9z 9£ 


on £ = 0 . 


Thus, the authors actually solved for a constant stress. However, this constant turns out to be zero for this 
solution. The boundary conditions and solutions for the E ^ 2 layer are given in Appendix A. 

Now that a characteristic velocity, U, has been determined, the values for e and M can be calculated 
and are given in Table 3. It is noted that e is not small. Since, for the linear results to be valid, e must be 
small, the question of the value of the linear work is raised. For the linear floating zone problem being 
studied, analytical solutions can be found and this means that there is a complete description of the flow 
dynamics. In general, analytical solutions cannot be found for strongly nonlinear problems. The great value 
of linear results is that they provide a frame of reference for analyzing the nonlinear contributions. Nonlinear 
effects can often be seen as secondary influences on the basic linear flow even for relatively large values of 
the parameter estimating the nonlinearity. The linear results were found to be valuable for interpreting the 
weakly and strongly nonlinear results presented in this paper. Also, the results presented in this paper reveal 
that e may not be so large for a more realistic model, thus reducing the nonlinear effects. 


3.5 The Ekman Layers Associated with the E */ 2 Layer 

Tire E^/3 layer azimuthal flow, Vj [equation (A-13)], does not vanish on the endwalls. Ekman 
layers are required on these boundaries to satisfy the no-slip condition. This report now proceeds to examine 
these layers. A transformation to f = E ' 1 / 2 (1-z) is made. This discussion is limited to the Ekman layer at 

f = 1. [A similar discussion for the Ekman layer at z = -1 is pursued using the transformation, £' = E ‘^/ 2 
(1 + z).] Taking v i /2 ~ ^0) to match Vj^ and satisfy the no-slip condition on the endwall, the following 
amplitude scales are obtained for the Ekman layer: 
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u 1/2~0(D , w 1/2 ~e'/ 6 , P 1/2 ~E 2 « , ^1/2-E 1 / 2 . 

Since Uj^ > u l/3 the appropriate boundary condition is that vl\J2 ~ 0 on f = 0. The equations, boundary 
conditions, and solutions for this layer are given in Appendix B. 

3.6 The First Order E l/ 2 Sidewall Boundary Layer and the Sidewall Boundary Layer 

This section now proceeds systematically to improve the accuracy of the results. The largest 
neglected term is in the boundary condition, w = 0, on the endwalls. One has w i /3 = 0, but w i /2 ~ ^ • 

Therefore, a first order E^' 2 sidewall boundary layer, with all the previous scales multiplied by the factor 
E^/6, i s introduced. These amplitude scales are: 

l u l/3 ~ E 1/2 , JVJ/3-E 1 / 6 , JWJ/3-E 1 / 6 , !P 1/3 -H 1 / 2 , ^ 1/3 ~E 1 / 2 . 

The preceding suffix 1 denotes a first order quantity. The equations are given in Appendix C. 

Ekman compatibility yields, 

1 ^ 1/3 = ±Ivi/ 3 , on z = ±l . (19) 

To proceed there is no need to solve for the first order E^ 2 layer, but the azimuthal stress boundary 

condition must be examined. Since the leading order E^ 2 layer satisfies the vanishing azimuthal stress condi- 
tion, an attempt is made to apply the same condition to the first order layer. Consider the following inte- 
grated form of equation (C-7), 


+ 1 2 ° 9 2 


a 


d£ dz = 2 


-1 0 


d£ dz 


f + ^ 2 1^1/3 r a ~ ~ i 

"J i f=0 dZ = 2 J 0 l z =l “ 1^1/3 l z= _i)di . 

Using equation (19), one obtains, 



( 20 ) 
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The left-hand side of equation (20) is the integrated azimuthal sidewall stress and the right-hand side is not 

zero, hence, the first order E^ 3 i a y er cannot satisfy a vanishing azimuthal stress condition. A further 
boundary layer is required. 

Sidewall layers of thickness E^ 4 are well known in the rotating fluids literature and are usually 

invoked to satisfy azimuthal flow boundary conditions [9], An E 1 / 4 layer is introduced to satisfy the 

azimuthal stress requirement. It is known from the previous work on E^ 4 layers that the sidewall azimuthal 

stress is constant [see also equation (D-l 1)], independent of z. Thus, the first order E^ 3 layer azimuthal 
stress is taken to be independent of z, i.e., 


d 

dl 



= 0 


on f = 0 , 


which can be written as 


a 5 i^i /3 


on £ = 0 . 


( 21 ) 


This discussion on the formulation of the first order E 1 / 3 layer problem is completed by remember- 
ing that the leading order sidewall layer satisfies the no radial flow boundary condition and the axial 
thermocapillary stress condition, therefore, 

a2 l^l/3 

1 ^1/3 = 5 — = 0 » on £ = 0 , 

d£ 2 

is required. The complete set of boundary conditions is given in Appendix C. 

To formulate the E^ 4 sidewall layer equations and boundary conditions, one transforms to a 
boundary layer coordinate, 17, using t? = E _1 / 4 (A-r). The amplitude scales of the E 1 / 4 layer are determined 
by equating the magnitudes of the first order E 1 / 3 layer and E 1 / 4 layer azimuthal stresses. One obtains 

uj/q-E 7 / 12 , v 1 / 4 ~E 1 / 12 , w 1 / 4 ~E 1 / 3 , p j /4 ~ E 1 / 3 , t// 1/4 ~E 7 / 12 . 

This report now proceeds to find the expression for the first order E 1 / 3 layer azimuthal stress and 
hence for the E 1 / 4 layer azimuthal stress. The required boundary condition is, 

3 l v l/3 3 ^l/4 

— r— + — — , on £ = 77 = 0 . (22) 
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From equation (A-l 2) one obtains, 


oo oo 

/ vi/3L-,«=Zh)-' 

0 n=l 



( 23 ) 


Substituting equation (23) into equation (20) one obtains, 


r 3 i 7 i/3 


dz 


£=o 


: - 2 £(-D n+1 

n=l 



(24) 


Using equation (21), equation (24) can be integrated to give. 


9 l^l/3 

a? 


OO 

=-£(-i ) n+1 

^=0 S 



(25) 


Finally, from equations (22) and (25) one obtains, 


97 l/4 

3i? 


= f (-D n+1 

t?=° n =l 



(26) 


Note that the azimuthal stress boundary condition has been found without solving for the first order 
E 1 / 3 layer. 

The equations, boundary conditions, and solutions for the E^ layer are given in Appendix D. The 

boundary conditions and solutions for the Ekman layer associated with the E^ layer are also given in 
Appendix D. 

An attempt was made to solve for the first order E^ 3 sidewall layer and solutions involving multiple 
summation were obtained. The convergence properties of these forms presented substantial difficulties to 
the numerical computation of the solutions. Also, since accurate solutions for the linear problem were 
obtained using numerical methods (see Section 4) and since the essential dynamics of the linear flow are 

revealed without finding complete solutions (see Section 6), no analytical solutions for the first order E^ 3 
sidewall layer are given. 
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3.7 The Composite Analytical Solutions 


The dimensional analytical results for the leading order E 1 / 3 layer and its associated Ekman layer 

with the addition of the E 1 / 4 layer and its associated Ekman layer are shown in Figures 3. Figures 3 include 
the stream function, azimuthal velocity, temperature, axial velocity, radial velocity, and pressure. It should 

be noted that these plots do not include the first order E^ 3 layer. These results are discussed further in 
Sections 5 and 6. 


AXIS 

i 

AXIS I 




(a) Stream Function (cm 2 /sec); 1.10 x 10' 4 ; 
-1.10 x 10' 4 ; 2 x 10' 5 . 


(b) Azimuthal velocity (cm/sec); 3.51 x 10‘ 4 ; 
-2.76 x 1 O' 3 ; 3 x 10' 4 . 


Figure 3. The dimensional analytical solutions for the dependent quantities. Each solution consists 
of the sum of the leading order E 1 / 3 layer and its associated Ekman layer and the E 1 / 4 layer and 

its associated Ekman Layer. The first order E 1 / 3 layer is not included. To facilitate comparison 
with the linearized numerical results (Figs. 4), the analytical results are calculated for the 

same parameter values (AT = 5 x 10^°C, A = 1, E = 0.35 x 10' 2 , Pr = 0.023). The 
dependent variables are listed below each figure together with the dimensional 
maximum and minimum values and the contour increment. For the pressure, 
only the contour increment is given. 
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IV. NUMERICAL SOLUTIONS 


4.1 The Numerical Model 

In Section III it was found that the dimensionless parameter, e, multiplying the nonlinear terms in 
the governing equations, is not small for the parameters and model chosen. Thus, the linearized results, 
although valuable for understanding the flows and the flow confinement dynamics (see Subsection 3.4), 
cannot be considered as accurate. A numerical model was used to include nonlinear effects in these studies. 

For this work a code already developed for the Atmosphere General Circulation Experiment (AGCE) 
Program [11] was modified. The code is discussed in some detail by Fowlis and Roberts [12] and additional 
applications of the code, including the floating zone application, are described by Roberts et al. [13]. The 
numerical methods used in the code are described briefly in Appendix E. The original AGCE code was 
developed for a spherical shell configuration with inner and outer spherical boundaries and radial, latitudinal 
boundaries. The conversion to cylindrical geometry was carried out by making the radius of the inner sphere 
very large and considering only the small latitudinal range from the axis of symmetry to the outer latitudinal 
boundary which becomes the floating zone sidewall. The boundary conditions on the sidewall were modified 
to allow for a free surface in the presence of the axial thermocapillary stress. 

The code was thoroughly validated in its cylindrical form using accurate and detailed experimental 
measurements obtained using laser-Doppler velocimetry. The validations studied used the experimental data 
presented by Wam-Vamas et al. [14] and Hyun et al. [15,16] and excellent agreement was obtained for 
both linear and nonlinear flows. 


4.2 Numerical Results 

All the numerical results presented in this subsection were examined to be certain that the spatial 
resolution was satisfactory and that convergence occurred. Initially, a nonuniform grid of 25 (radial mesh 
points) x 40 (axial mesh points) was used, but, if the results indicated that more resolution was required, 
the grid resolution was systematically increased. The maximum number of grid points used was 40 x 120 
and with this resolution even the strongly nonlinear, unsteady flow was resolved. The iterative method for 
rapid convergence to a steady-state was used (see Appendix E). However, to be certain that a steady-state 
was the correct solution, the results were confirmed using uniform time-stepping. The unsteady flows were 
determined using time-stepping. As a consequence of the above and the validation studies described in 
Subsection 4.1, the numerical solutions presented are considered to be accurate solutions of the governing 
two-dimensional equations [equations (1) to (5)]. However, for nonlinear systems the solution computed is 
not necessarily unique. Different initial conditions can lead to different solutions (hysteresis effects). 
Further, as the nonlinearity is increased the flow is more likely to become unstable to azimuthal perturba- 
tions. 


To check the analytical results, an effectively linear flow was first computed. This was done by 

reducing AT from 5°C to 5 x 10~^°C, thus reducing e to 5.8 x 10"^, and keeping the same values for the 
other dimensionless parameters (see Table 3). The dimensional results for stream function, azimuthal 
velocity, temperature, axial velocity, radial velocity, and pressure are shown in Figures 4. A nonuniform grid 
of 25 x 40 and the iterative method for convergence were used. The linear flow results are discussed further 
in Sections V and VI. 

To examine the effects of nonlinearity, the code for AT = 0.5°C, corresponding to e = 5.8, was run, 
again keeping the same values for the other dimensionless parameters (see Table 3). The dimensional results 
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AXIS 



(a) Stream Function (cm^/sec); 3.74 x lCT^ 
;-3.74 x 10‘ 5 ; 1 x 10' 5 . 


(b) Azimuthal Velocity (cm/sec); 3.56 x 10" 4 ; 
-2.26 x 1 0" 3 ; 3 x 10' 4 . 


Figure 4. The dimensional solutions for the linearized numerical problem. This problem was linearized 
by taking AT = 5 x 10’ 4 °C and hence, e = 0.58 x 10'^. The other parameters have the values, 

'y 

A = 1, E = 0.35 x 10 , Pr = 0.023. The dependent variables are listed below each figure together 
with the dimensional maximum and minimum values and the contour increment. For the 
pressure, only the contour increment is given. 
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for stream function, azimuthal velocity, temperature, axial velocity, radial velocity, and pressure are shown 
in Figures 5. The results presented were computed with a nonuniform grid of 40 x 120 and time-stepping for 
about 60 sec from an initial state of no motion relative to the rotating frame of reference. However, the 
iterative method gave identical results. Note that the flow patterns are significantly different from the 
linearized results of Figures 4. These results are discussed further in Section 6. 


Finally, the code for AT = 5°C, corresponding to e = 58, was run, keeping the same values for the 
other parameters (see Table 3). This strongly nonlinear flow did not achieve a steady-state but settled to a 
flow with a weak oscillation. The dimensional results for two extreme phases of the oscillation are shown 
in Figures 6 and 7; stream function, azimuthal flow, and temperature plots are presented. These results were 
computed with a nonuniform grid of 40 x 1 20 and time-stepping for about 40 sec from an initial state of no 
relative motion. This flow was also examined starting from the same initial state but using the iterative 
method for a time estimated as 20 sec and then time-stepping for another 25 sec. The solutions achieved 
were essentially the same as those presented in Figures 6 and 7. These results are discussed further in Section 


VI. 

AXIS 

i 



(a) Stream Function (cm 2 /sec); 8.10 x 10 2 ; 
-8.10 x 10‘ 2 ; 1.5 x 10' 2 . 


AXIS 



(b) Azimuthal Velocity (cm/sec); 0.351 ; 
-0.417; 0.1. 


Figure 5. The dimensional solutions for the weakly nonlinear flow. For this problem AT = 0.5°C 

and hence, e = 5.8. The other parameters have the values, A = 1 , E = 0.35 x 10 , Pr = 0.023. 
The dependent variables are listed below each figure together with the dimensional maximum 
and minimum values and the contour increment. For the pressure, only the contour 

increment is given. 
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AXfS 



(a) Stream Function (cm^/sec); 0.846; -0.840; 0.15. 


(b) Azimuthal Velocity (cm/sec); 2.65; -0.432; 0.5. 


AXIS 

I 



(c) Temperature (°C); 5.00; 0; 0.5. 

Figure 6. The dimensional solutions for the strongly nonlinear flow. For this problem AT = 5°C 

and hence, e = 58. The other parameters have the values, A = 1, E = 0.35 x 10"^, Pr = 0.023. 
This figure shows an extreme phase of the observed oscillation; Figure 7 shows the other 
extreme phase. The dependent variables are listed below each figure together with the 
dimensional maximum and minimum values and the contour increment. 
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AXIS 


AXIS 




(a) Stream Function (cm 2 /sec); 0.840; -0.847; 0.15. (b) Azimuthal Velocity (cm/sec); 2.66; -0.435; 0.5. 

AXIS 

I 



Figure 7. The dimensional solutions for the strongly nonlinear flow. For this problem AT - 5°C 

-9 

and hence, e = 58. The other parameters have the values, A = 1, E = 0.35 x 10 , Pr = 0.023. 
This Figure shows an extreme phase of the observed oscillation; Figure 6 shows the other 
extreme phase. The dependent variables are listed below each figure together with the 
dimensional maximum and minimum values and the contour increment. 
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V. COMPARISON OF THE ANALYTICAL AND LINEARIZED NUMERICAL RESULTS 


The analytical results (Fig. 3) are now compared with the linearized numerical results (Fig. 4). It 

should be remembered that Figures 3 do not include the first order E 1 / 3 layer and that some higher order 
flows do not satisfy the boundary conditions. 

Very good agreement between the analytical and numerical results is obtained for the azimuthal 
flow in both overall contour shapes and maximum and minimum values (see Figs. 3b and 4b). For the 

analytical results the azimuthal stress due to the E*/ 3 layer vanishes on the sidewall. The combined 
azimuthal stress due to the first order E^ 3 layer and the E^/ 4 i a y er a i so vanishes on the sidewall but, of 
course, the first order E 1 / 3 layer is not included in Figure 4b. The amplitude of v 1/3 is 0(1) and the omitted 
v l/3 ~ (0.39) is less than vj ~ E^ 2 (0.62). 

The numbers in parentheses after the amplitude scales show the magnitudes of the amplitude scales 

for the particular variables. These values are based on E = 0.35 x 10’ 2 (see Table 3). Since E is not exceed- 
ingiy small, the amplitudes just discussed, and those in the discussion to follow, are sometimes not very 
different. Thus, these amplitude comparisons should not be taken as truly quantitative, since factors in the 
exact solutions could make significant contributions. 

Good agreement is obtained for the axial flow in both overall contour shape and maximum and 
minimum values (see Figs. 3d and 4d). For the analytical results, the E 1 / 3 layer satisfies the axial thermo- 
capillary stress. However, the axial stress due to the E 1 / 4 layer does not vanish on the sidewall. The ampli- 
tude °f w ] is 0(1) but the omitted j Wj^ 3 ~ E^ (0.39) is not less than Wj/ 4 ~ E^/ 3 (0.15). 

Good agreement is obtained for the radial flow in overall contour shape but there is not good agree- 
ment for flow magnitudes (see Figs. 3e and 4e). Note that the extreme negative values for the radial flow 
occur in the sidewall comer regions and that these regions are not resolved in the analytical work. For the 

analytical results the radial flow due to the E^ 3 layer vanishes on the sidewall [see equation (A-13)] but 
the radial flow due to the E 1 / 3 layer, Ekman layer does not vanish on the sidewall [see Fig. 4e and equa- 
tion (B-10)] . Also, the radial flow due to the E 1 / 4 layer does not vanish on the sidewall [see Fig. 4d and 
equation (D-10)]. The amplitude of u }/3 is OfE 1 / 3 ) (0.15) but the omitted jUj^ ~ E 1 / 2 (0.059) is not 
less than uj/ 4 ~ E 7 / 12 (0.037). Good agreement is obtained for the pressure (see Figs. 3f and 4f). The 
amplitude of p j y 3 is 0(E 1 / 3 ) (0. 1 5 ) and the omitted j P i / 3 ~ E 1 / 2 (0.05 9) is less than p j ~ E 1 / 3 (0. 1 5 ). 
Poor agreement is obtained for the stream function (see Figs. 3a and 4a). For the analytical results this is a 
consequence of the radial flow at the sidewall due to the E 1 / 4 layer. The amplitude of i//^ 3 is OtE 1 / 3 ) 
(0.15) but the omitted ^]/ 3 ~ E 1 / 4 (0.059) is not less than ^jy 4 ~ E 7 / 12 (0.037). To reduce the E 1 / 4 

layer radial flow at the sidewall to zero, it is necessary to introduce a second order E 1 / 3 layer. The analytical 

studies were not pursued beyond this point. As expected, the temperature results are identical (see Figs 3c 
and 4c). 
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VI. DISCUSSION OF RESULTS 


6. 1 The Linear Flow 

The agreement between the analytical and linearized numerical results and the authors understanding 
of the discrepancies between these results (see Section V) are such that it is believed the analytical results 
reveal the essential dynamics of the linear flow. Since the linear flow is symmetrical about the mid-depth 
plane of the floating zone, the following discussion is confined to one half of the full floating zone. The 
surface tension of liquid silicon increases with decreasing temperature, hence, the sidewall stress pulls fluid 
away from the mid-depth plane towards the endwalls. In the vicinity of the endwall most of this fluid turns 

around before entering the Ekman layer and circulates in an E^^ sidewall layer. As a consequence of the 
basic rotation, the inward radial flow near the endwall leads to positive azimuthal flow and the outward 
radial flow near the mid-depth plane leads to negative azimuthal flow. The azimuthal flow near the endwall 

drives an Ekman layer on that wall over the E layer. The Ekman layer has a meridional flow which 

circulates through a first order E^/^ layer. The flow from the first order layer enters the Ekman layer in 
the comer region, not resolved by the analytical studies, then returns directly from the Ekman layer. 

The leading order E^ layer satisfies all the sidewall boundary conditions, but the first order layer 
cannot satisfy the sidewall azimuthal stress condition. Thus, an E^ layer with negative azimuthal flow 
exists in which the meridional flow circulates in the opposite sense. In combination with the first order E^ 
layer, this E^ i a y er satisfies the vanishing azimuthal stress condition. The high order radial flow and the 
axial stress introduced by the E^ layer remain unsatisfied and require the presence of a still higher order 
E^/3 i a yer (or layers). 

In agreement with what was found by Smith [7] and Smith and Greenspan [8], these linearized 
results (see Fig. 4) show, for the weak thermocapillary flow, that rotation does confine the flow to sidewall 
boundary layers and that there is no flow in the interior of the melt relative to the rotating endwalls. For 
comparison, the flow for identical conditions to the flow of Figures 4, except that the rotation rate was 
made zero, is shown in Figures 8. Figures 8 show the stream function, the axial velocity, and the radial 
velocity. The azimuthal flow is not shown since it is zero for a nonrotating flow and the temperature is not 
shown since it is virtually identical to the conduction solution (see Fig. 4c). [The nonrotating results were 
computed for a nonuniform 40 x 70 mesh and the iterative method (see Appendix E) was used.] It is seen 
that this flow is considerably stronger than the rotating flow and that it penetrates much more deeply into 
the interior of the melt. 

For the half zone. Smith [7] and Smith and Greenspan [8] did not require an E^ ] a yer. Thus, the 
presence of this thicker layer in the full zone means that confinement is not quite as good as for the half 
zone. The boundary layers could be made thinner by increasing the rotation rate, £2, and hence decreasing 
E. 


Increasing will also increase the flow due to centrifugal buoyancy and it should be remembered 
that this flow has been neglected in the present study. Smith [7] and Smith and Greenspan [8] studied 
the linearized flow due to centrifugal buoyancy in the half zone. They found a flow which circulated 
through the interior via Ekman layers on the endwalls. But, for the conditions given in Tables 1 and 2, they 
found this flow to be much weaker than the thermocapillary flow. For the full zone this flow will be still 
weaker, since in this configuration one cannot have Ekman layers at the mid-depth plane allowing for radial 
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AXIS 



AXIS 



(a) Stream Function (cm 3 /sec); 3.24 x 10' 4 ; 
-3.24 x 10 A ;S x 10‘ 5 . 


(b) Axial Velocity (cm/sec); 4.04 x 10‘ 3 ; 
-4.04 x 1 0‘ 3 ; 5 x 10* 4 . 


AXIS 



(c) Radial Velocity (cm/sec); 1.29 x 10‘ 3 ; 
-1.43 x 1 O' 3 ; 3 x 10 -4 . 


Figure 8. The dimensional solutions for the nonrotating flow and for AT = 5x 10" 4 °C The 
other parameters have the values, A = 1, Pr = 0.023. The dependent variables are listed below each 
figure together with the dimensional maximum and minimum values and the contour increment. 


24 



fluid motion. Another parameter to be considered as F2 is increased in the rotational Bond number, B (see 
Section II). For the cylindrical geometry of this model to remain accurate B must remam small. 


For this work realistic values were chosen for a floating zone system. As a consequence, although E 
is small, it is not very small and the ratios between the amplitudes of the higher order and lower order 
boundary layer quantities are not particularly small, It could be argued that to check further this analytical 
theory, one should perform numerical studies for smaller values of E and look for systematic converge ot 
the analytical and numerical results. Some such studies were performed and convergence was observed. 
However, since the linearized flow is somewhat unrealistic and was performed for the purpose of under- 
standing more deeply the nonlinear flows, the authors chose not to belabor this point and hence the linear 
results for the more realistic value of E only are presented. 


6.2 The Nonlinear Flow 

The weakly nonlinear numerical results for e = 5.8 (AT = 0.5°C) are shown in Figures 5. The flow is 
now, of course, much stronger and substantially different flow patterns are revealed when compared with the 
linearized results of Figures 4. The meridional flow is no longer a simple sidewall boundary layer flow, but 
takes the form of two vortices displaced towards the endwalls. The stronger meridional flow also penetrates 
further into the interior especially near the endwalls. The azimuthal flow is negative everywhere on the 
sidewall and has positive maxima in the interior. The isotherms reveal some slight distortion from the con- 
duction solution (Fig. 3c) informing us that nonlinear temperature advection is becoming significant. 
Although the influence of rotation in confining this flow is clearly reduced, some confinement is still 
present. 


Solutions of the strongly nonlinear flow for e = 58 (AT = 5°C) are shown in Figures 6 and 7. The 
flow is rapid and unsteady. A weak oscillation is present and Figures 6 and 7 show extreme phases of the 
oscillation. The instantaneous flow is not symmetrical about mid-depth, but the time average is symmetrical. 
The oscillation period is about 0.4 sec. The meridional and azimuthal flows penetrate deeply into the 
interior. The maximum axial velocity fluctuates around 12.5 cm/sec and occurs on the sidewall. The 
maximum radial flow is inward; it fluctuates around 7.75 cm/sec and occurs near the comer regions. The 
conduction temperature solution (Fig. 3c) is severely distorted by the strong flow. 

This flow is very nonlinear and the comments made in Subsection 4.2 about uniqueness are relevant. 
The flow presented in Figures 6 and 7 may not be unique and it is almost certainly unstable to azimuthal 
perturbations. Time dependent flows are believed to be the cause of dopant striations in crystals (see Section 
I) and theoretical and experimental studies of the stability of thermocapillary flows have been performed 
[17,18,19,20]. 


VII. CONCLUSIONS 


For the model and parameter values chosen, the present work establishes that the thermocapillary 
flow is strong and that rotation cannot confine the flow. This conclusion is apparently reinforced by the 
fact that AT = 5°C could be considered too small and hence real flows are still stronger. However, an assess- 
ment of the results and the model reveals that actual flows may not be so strong and that rotation could 
still be used for confinement. 

The distortion of the isotherms for the strongly nonlinear flow shown in Figures 6 and 7 reveals 
that temperature advection by the flow is significant even for the small Prandtl number (Pr = 0.023) of 
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silicon. In the authors model the temperature on the sidewall was specified but, in fact, the actual tempera- 
ture on the sidewall is determined by a balance of radiation, temperature advection, and conduction and 
since the results in this report show that temperature advection is important, a more realistic temperature 
boundary condition will be obtained by including all the processes. The inclusion of temperature advection 
will reduce the temperature gradient on the sidewall away from the endwalls and increase the temperature 
gradient close to the endwalls. This conclusion was also reached by Sen and Davis 119] who examined 
thermocapillary flow in a simpler configuration. Since the thermocapillary driving stress is now weaker 
over the bulk of the fluid and since viscosity may reduce the flow close to the endwalls, the overall non- 
linear flow may be substantially reduced. A reduction in the flow strength will permit the confining effect 
of rotation to play a stronger role. 


Since the melting point of silicon is about 141 0°C, it is not unreasonable to think of AT = 5”C as 
somewhat small. However, since crystals have sharp melting points, the whole floating zone system could 
be placed in a furnace with a background temperature just below the melting point and only a small amount 
of heat from the ring heater would then be required to produce a floating zone [22, 23]. This would result 
in a small AT and a corresponding weak flow. Surfactants will also reduce the surface tension stress [24,25] 
but, of course, suitable surfactants, which could not be considered as contaminants, would have to be used. 

It is concluded that rotation may still be a means of confining the thermocapillary flow in a floating 
zone system to the melt sidewall. If this is so, then, in space, all flows in the melt interior will be greatly 
reduced. Further theoretical work and experimental confirmation of the theory is needed before this sugges- 
tion can be answered conclusively. 
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APPENDIX A 


THE LEADING ORDER E 1/3 SIDEWALL BOUNDARY LAYER 


A transformation to the boundary layer scale was made using, 


£ = E -1 / 3 (A-r) . 


The velocity scale U was chosen to make Wj / 3 of order unity. The amplitude scales are: 
u l/3 = e1/3 u 1/3 ’ v l/3 = ^1/3 > w l/3 = ™l/3 , P 1/3 = E 1/3 P 1/3 , ^ 1/3 


where the tilda denotes an 0(1) quantity. The scaled governing equations are, 


2 v 


1/3 


a Pl/3 

3£ 


9 ^1/3 _ 

— ~2u l/3 = 0 , 

a 2 wi / 3 3p 1/3 




9z 


Suj/3 a ^l/3 

9£ + ~dT~ ~ ° ’ 

In terms of the stream function, one has, 


l l/3 


9 ^ 1/3 


, w, 


3^1/3 


9z ’ 1/ 3 3£ 

The eliminated equation for $ 1/3 is. 


9 6 ^l/3 3 2 ^ 1/3 


+ 4- 


= 0 . 


9£ f 


9z^ 


The boundary conditions on 4>\j3 are: 


(A-l) 


^1/3 > 


(A-2) 


(A-3) 


(A-4) 


(A-5) 


(A- 6 ) 


(A-7) 
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(A-8) 


^ i /3 0 , as £ 

%/3 


d 2,f ' - 




^1/3 


+ sin Vi 7rz = 


1/3 




3$* 


= 0 , on £ = 0 , 


\[/ 1/3 = 0 , or z = ±l 


(A-9,a,b,c) 


(A- 10) 


1 he solutions are found using separation of variables. Sin Vi ttz in equation (A-9b) is Fourier analyzed 


as, 


sin Vi it z = aj, sin n 7r z , 
n=l 


(All) 


where 


a n = (-D n+1 — f 

7r(n 2-1/4) 


Note that equation (A-l 1) is satisfied for all z in the domain -1 < z < 1 except z - ±1 . 
The solutions are: 


~ v' a n T ~7nl 2 -Viy n ^ . / 3 27tY| . 

n=l z ^n 


sin mrz 


? i o - E- ^ [^" 5 + 1 e ' ,/27ni s ” (| ^)] 

n=l n 


cos nwz , 


w 


1/3 


e ~Tn^ 2 - 1/2 T, 
6 "3 e 


n?sin (|'?J)] si 


sin n7rz , 


oo *- 

n=l 11 


V" a n 7n T ”Tn^ , 2 - 1 / 2 T n ^ • / 3 fc 2ir \l 
“1/3 = 2--— [ e + 3 e “ sm (2 T n { 'T)J 


COS n7TZ 


n=l 


(A-l 2) 


(A-13) 


(A- 14) 


(A- 15) 


Pl/3 =X — \ [ e_Tn ^ /27n ^ Sin (l ~ y)] C ° S n7rZ ’ ( A ' 1 6) 

n=l T n 

where y n = (2n7r)^ Note that all the boundary conditions, equations (A-8) through (A-10), ^and the 
symmetry conditions, equations (12) and (13) are satisfied, and also that the azimuthal stress, Sv^/di;, 
vanishes on £ = 0. 
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APPENDIX B 


THE EKMAN LAYERS ASSOCIATED WITH THE E 1/3 LAYER 

A transformation to the boundary layer scale was made using, 

f = E~ 1 / 2 (l-z) . (B-l) 

Here only the Ekman layer near z = 1 is discussed. The Ekman layer near z = -1 is found by using ?' = 
1/9 ° 

E” ' ( 1 + z). Using the Ekman layer to satisfy vj / 3 + vj /2 = 0 at z = 1 , the amplitude scales are determined 

as: 


u l/2 “ “1/2 ’ v l/2“' r i/2 > w l/2 = E 1/6 wj/2 , Pj/ 2 = E 2/3 P 1/2 , *j/ 2 = E 1/2 * 1/2 • 
The scaled governing equations are: 


d 2 “l/2 

3f 2 


+ 2 ^ 1/2 


= 0 , 


a2 ^l/2 

9? 2 



= 0 , 


3 2 w ]/ 2 _ 3p 1/2 

a? 2 = "“ 


3u 1/2 ^3w 1/2 

3£ + 


(B-2) 


(B-3) 


(B-4) 


(B-5) 


The eliminated equation for u ^2 or ^l/2 


d4tr i/2 

ar 4 


+ 4ui/2 = ° . 


(B-6) 
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The boundary conditions are: 


All dependent quantities -*• 0 , as f -*■ 00 , 

(B-7) 

uj 12 = 0 , on f = 0 , 

(B-8) 

' r i/2 + h/3 = 0 * on $" = 0 • 

(B-9) 

Tlie solutions obtained are: 


= V(|)e'f (cosf + sin?) , 

(B-1U) 

vj / 2 = v (?) e ‘^ cos ? , 

(B-l 1) 

w 1 / 2 = l - V'«) e* (cos f + sin f ) , 

(B-12) 

“ 1/2 = v (£) e '^ sin ? » 

(B-l 3) 

Pi /2 = v '«) e * ?sin r . 

(B-l 4) 


where 


v«)=2>, n -^- 





and V'(£) is its derivative with respect to £. Note that the boundary conditions [equations (B-7) through 
(B-9)] are satisfied. Note also that and ™l/2 are nonzero on f = 0. The correct boundary conditions 

for these quantities are satisfied in conjunction with the first order layer. Note that Uj/ 2 i s no t equal 
to zero on £ = 0. This point is discussed further in Section 6. 
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APPENDIX C 


THE FIRST ORDER E 1 / 3 SIDEWALL BOUNDARY LAYER 


Recalling the amplitude scales for the leading order problem and our boundary layer requirement 
that, Wj /3 ~ El/**, th e dependent quantities are expanded as follows: 


u l/3 = E 1 / 3 ( u i/3 + E 1 / 6 iUj /3 + ... ) , (C-l) 

v l/3 = v l/3 + E 1 / 6 l v l /3 + — . (C 2) 

w 1/3 = wi/3 + E 1/6 1 w 1 / 3 + ... , (C . 3) 

Pl/3 = E 1/3 CPl/3 + E 1 / 6 1 P 1/3 + ... ) , (C-4) 

^1/3 = e1 / 3 (0j/3 + E 1 / 6 J0J/3 + ... ) . (C-5) 


Using this sidewall layer to match the Ekman layer w ^j2 ~ E^, one obtains for the amplitude scales. 
1 U 1/3 = E 1/2 1U1/3 , l v l/ 3 = E 1/6 !V 1/3 , iw 1/3 = E 1 / 6 jw 1/3 , 1P1/3 = E 1 / 2 jp 1/3 , 

1^1/3 = E 1/2 1^1/3 • 

The scaled governing equations are identical to those of the leading order problem, namely, 

- ~ _ d lPl/3 

2 \ v \I3 > (C-6) 

a Vl /3 

— ~ 2 2 1 u 1/3 “ 0 > ( C - 7 ) 

a2 l^l/3 _ djPj/3 

3^2 3z ’ (C-8) 
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(C-9) 


a,u 1/3 [ 9 i™i /3 
9£ 9z 

The stream function is written as, 

9 l^l/3 „ _ 9 l^l/3 

1 U 1 /3 9 z ’ l w l/3 ’ 


The eliminated equation for 1 ^ 1/3 is. 


9 6 1^1/3 a2 l^l/3 

+ 4 

9| 6 9z 2 


The boundary conditions for j \p j / 3 are: 


1 ^ 1 / 3^0 , as £ 


3 1^1/3 9 Vl/3 

l^]/3 = T- = 0 - on ?-0 , 

1 9£ 2 9^ 5 


l^l /3 = ± Ivi /3 = ± i v (0 , on z = ±l , 


where 

V(£) = £ (-1 ) n 2 ^- + J e 1/27n? sin (| T n^)] 

n=l n 

as before. Solutions of the above equations and boundary conditions are not given. 


(C-1U) 

(C-ll) 

(C-12) 
(C- 13) 
(C-14) 
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APPENDIX D 


THE E 1 / 4 SIDEWALL BOUNDARY LAYER AND ITS ASSOCIATED EKMAN LAYER 


A transformation to the boundary layer scale was made using, 
r? = E'l/ 4 (A-r) . 


To take care of the sidewall azimuthal stress introduced by the first order E 1 / 3 sidewall layer, 
OfE 1 / 12 ) is required and the amplitude scales are then determined as. 


u l/4 - E?/12 “1/4 > v l/4 = E ‘ /12 ' r i/4 ■ w J/4 = E 1 / 3 wj/4 , p 1/4 = E*/ 3 p, /4 , 

*l/4 = E 7 /> 2 ? 1/4 . 


The scaled governing equations are: 


2 vi / 4 = - 


3pi/4 

drj ’ 


d z v 


1/4 


hr} 1 


-2ui/ 4 = 0 , 


aF,/4 


3z 


= 0 , 


9UJ/4 3w 1/4 

•+ — = 0 




9z 


The stream function is written as, 


9^1/4 _ a^i /4 

U l/4 = . w l/4 = 


3r? 


(D-l) 
v l/4 ~ 


(D-2) 

(D-3) 

(D-4) 

(D-5) 

(D-6) 
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From equations (D-2) through (D-6) one deduces, 


auj /4 3vj/ 4 apj /4 3 2 w!/ 4 d 2 $\/4 

dz 3$ dz dz 2 dz 2 


(D-7a,b,c,d,e) 


One clearly has an Ekman layer associated with this E ^ layer. Scaling for this Ekman layer to make 
v l /4 + v l /2 = 0 on z = ±1, yields: 


u 1/2 = E 1/12 ui/ 2 , v 1/2 = E 1 / 12 7 1/2 , w 1/2 = E 1 / 3 w 1/2 , p = E 5 / 6 p 1/2 , * 1/2 = E V/12 5 1/2 • 


Ekman compatibility yields. 


w 



b l±!± 

dr) 


on z = ±1 


(D-S) 


Equation (26) gives the azimuthal sidewall stress and let this stress equal S, i.e., 


97 l/4 

3t? 


=£ 

n=l ?n 


s 


on rj = 0 


(D-y) 


Remembering that the problem is symmetric about z = 0 and using 

the above equations and 

boundary conditions, the solutions are found to be. 

^l/4 = -^ Se ' T?z > 

(D-10) 

v 1/4 = -Se- T ? , 

(D* 11) 

1 „ 

w l/4 = 2 Se ^ z > 

(D-12) 

ir i/4 = -^ Se ' T? > 

(D-13) 

p l/4 =-2Se~ T l . 

(D-14) 

Note that the azimuthal sidewall stress condition is satisfied but i/\j 4 and u ^4 

are not equal to zero on 


T? = 0. 
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The Ekman layer solutions are found 
order boundary conditions are: 

as before for the E^ layer (see Appendix B). The leading 

o 

ii 

(N 

1 

1 on f = 0 

(D-15) 

7 l/2 +V l/4 = 0 • , 


(D-16) 

The solutions are: 



«/> i /2 = j Se' 7 ? e'f (sin f + cos f ) , 

(D-17) 

v, j 2 ~ Se' 7 ? e'f cos f 

5 

(D-18) 

wj /2 = - j Se* 7 ? e'f (sin ? + cos f) , 

(D-19) 

uj j 2 = Se* 7 ? e'^ sin f , 


(D-20) 

Pj /2 = Se'^e'f sin f . 


(D-21) 


These Ekman layer solutions rigorously satisfy the E layer endwall boundary conditions. 
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APPENDIX E 


THE NUMERICAL CODE 


The code solves the Boussinesq form of the axisymmetric, incompressible, Navier-Stokes equations 
in the configuration of Figure 2. Finite difference methods are used. To give adequate resolution in the 
boundary layers without wasting resolution in the interior, a nonuniform computational mesh is specified. 
The dependent variables are represented by their values at the same primary mesh points, which simplifies 
part of the coding and allows for the efficient use of implicit methods. The spatial representation conserves 
several integral quantities. 

An iterative method based on time-stepping is used to obtain steady solutions. This method has the 
option of different time steps for temperature and velocity and these steps vary smoothly according to 
specified powers of the mesh spacing in the two directions. This permits more rapid convergence to steady 
states than simple uniform time-stepping. An alternating-direction implicit formulation is used for advection 
and diffusion and for the Coriolis force, with implicit parameters chosen to damp high frequency 
phenomena. 

The pressure field is determined from the continuity equation. The total pressure at the previous 
iteration is used as a first approximation to update the velocity variables. The velocity divergence is then 
calculated and used to obtain a convection pressure field which is then used to correct the velocities. This 
convection pressure satisfies a representation of a Poisson-like equation, which is solved by alternating direc- 
tion implicit iteration. 


37 





1. REPORT NO. 

NASA TP-2576 

2. GOVERNMENT ACCESSION NO. 

3. RECIPIENT'S CATALOG NO. 

4. title and subtitle 

Analytical and Numerical Studies of the Thermocapillary Flow 
in a Uniformly Rotating Floating Zone 

5. REPORT DATE 

March 1986 

6. PERFORMING ORGANIZATION CODE 

7. AUTHOR(S) 

William W. Fowlis and Glyn 0. Roberts* 

0. PERFORMING ORGANIZATION REPORT # 

9. PERFORMING ORGANIZATION NAME AND ADDRESS 

George C. Marshall Space Flight Center 
Marshall Space Flight Center, Alabama 35812 

10. WORK UNIT NO. 

M-518 

1 1. CONTRACT OR GRANT NO. 

13. TYPE OF REPORT & PERIOD COVERED 

Technical Paper 

14. SPONSORING AGENCY CODE 

12. SPONSORING AGENCY NAME AND ADDRESS 

National Aeronautics and Space Administration 
Washington, D.C. 20546 


15. SUPPLEMENTARY NOTES 

Prepared by Space Science Laboratory, Science Engineering Directorate. 


* Roberts Associates, Inc., Vienna, Virginia 22180. 

16. ABSTRACT _ ' ^ ™“ ~ ~ ~ ~ ~ “ 

The microgravity environment of an orbiting vehicle permits crystal growth experiments in the presence of greatly 
reduced buoyant convection in the liquid melt. Crystals grown in ground-based laboratories do not achieve their potential 
properties because of dopant variations caused by flow in the melt. The floating zone crystal growing system is widely used 
to produce crystals of silicon and other materials. However, in this system the temperature gradient on the free sidewall 
surface of the melt is the source of a thermocapillary flow which does not disappear in the low-gravity environment. 

Smith and Greenspan theoretically examined the idea of using a uniform rotation of the floating zone system to 
confine the thermocapillary flow to the melt sidewall leaving the interior of the melt passive. These workers considered a 
cylinder of fluid with an axial temperature gradient imposed on the cylindrical sidewall. They considered a half zone and 
examined the linearized, axisymmetric flow in the absence of crystal growth. They found that rotation does confine the 
linear thermocapillary flow. 

In this paper the simplified model of Smith and Greenspan is extended to a full zone and both linear and non- 
linear thermocapillary flows are studied theoretically. Analytical and numerical methods are used for the linear flows and 
numerical methods for the nonlinear flows. It was found that the linear flows in the full zone have more complicated and 
thicker boundary layer structures than in the half zone, and that these flows are also confined by the rotation. However, 
for the simplified model considered and for realistic values for silicon, the thermocapillary flow is not linear. The nonlinear 
flows were examined by first computing a weakly nonlinear flow and then computing the fully nonlinear flow. The weakly 
nonlinear flow is steady, has less boundary layer character, and penetrates more deeply into the interior than the linear 
flow but still shows some rotational confinement. The fully nonlinear flow is strong and unsteady (a weak oscillation is 
present) and it penetrates the interior. Some non-rotating flow results are also presented. 

Since silicon has a large value of thermal conductivity, one would expect the temperature fields to be determined 
by conduction alone. This is true for the linear and weakly nonlinear flows, but for the stronger nonlinear flow the results 
show that temperature advection is also important. Thus, this work reveals that for the nonlinear flow, a radiative sidewall 
boundary condition would be an improvement over the specified temperature boundary condition used in this paper and 
previously by others. Such a boundary condition would weaken the sidewall axial temperature gradient and hence the 
thermocapillary flow allowing the confining effect of rotation to play a stronger role. Hence, uniform rotation may still 
be a means of confining the flow and the results obtained define the procedure to be used to examine this hypothesis. 


V7. KEV WORDS | 18.. DISTRIBUTION STATEMENT 


Floating Zone Unclassified — Unlimited 

Theoretical Thermocapillary Flow Studies 
Confinement of Thermocapillary by Rotation 


19 . SECURITY CLASSIF. ( of thU report* 

Unclassified 


Subject Category 34 


20. SECURITY CL ASSI F. (of this page) 

21. NO. OF PAGES 

22. PRICE 

Unclassified 

41 

A03 


For sale by National Technical Information Service, Springfield, Virginia 22161 

NASA-Langley, 1986 




















National Aeronautics and 
Space Administration 
Code NIT-4 

Washington, D.C. 
20546-0001 


OHicial Business 

Penalty lor Private Use, S300 





POSTMASTER: 


If Undelhrerable (Section 158 
Postal Manual) Do Not Return 




